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AN  OPTIMAL  DOMAIN  DECOMPOSITION  PRECONDITIONER 

FOR  THE  FINITE  ELEMENT  SOLUTION 

OF  LINEAR  ELASTICITY  PROBLEMS 

BARRY  F.  SMITH* 

Abstract.  For  linear  elasticity  problems  the  finite  element  method  is  an  extremely  successful 
method  to  model  complicated  structures.  The  successful  implementation  requires  the  solution  of  very 
large,  sparse,  positive  definite  linear  systems  of  algebraic  equations.  A  new  technique  for  solving  these 
systems  using  preconditioned  conjugate  gradients  is  proposed.  Using  ideas  from  both  additive  Schwarz 
methods  and  iterative  substructuring  methods,  we  prove  that  the  condition  number  of  the  resulting 
system  does  not  grow  as  the  substructures  are  made  smaller  and  the  mesh  is  refined.  This  result  holds 
for  two  and  three  dimensions.  Numerical  experiments  have  been  performed  to  demonstrate  the  power 
of  this  method.  For  linear  elasticity  problems  in  two  dimensions  the  condition  numbers  are  observed 
numerically  to  be  less  than  four  when  using  a  regular  mesh. 
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1.  Introduction.  We  introduce  a  new  domain  decomposition  method  of  solving 
the  Hnear  systems  arising  from  the  finite  element  discretization  of  coupled  elliptic 
systems  using  a  preconditioned  conjugate  gradient  method.  The  method  incorpo- 
rates certain  features  of  both  additive  Schwarz  methods  and  iterative  substructuring 
methods.  We  prove  that  for  second  order,  self  adjoint,  uniformly  elliptic  systems  the 
condition  number  of  the  resulting  preconditioned  linear  system  is  bounded  indepen- 
dently of  the  number  of  substructures  and  the  mesh  refinement.  The  method  is  highly 
parallelizable. 

Very  early,  Sobolov  [31]  showed  that  the  Schwarz  alternating  method  converges 
for  the  equations  of  linear  elasticity.  More  recent  work  using  the  Neumann- Dirchlet 
algorithm  was  done  by  Bj0rstad  and  Hvidsten  [4].  De  Roeck  [15]  has  implemented 
an  iterative  substructuring  type  algorithm  for  elasticity  problems.  Hughes  and  oth- 
ers have  been  using  element-by-element  preconditioning  [22], [33],  on  large  structural 
problems.  The  preconditioned  problems  for  these  latter  two  algorithms  still  can  re- 
quire hundreds  of  conjugate  gradient  iterations. 

For  the  p-version  finite  elements  Mandel  has  analyzed  iterative  substructuring 
type  algorithms  for  elasticity  [25], [26], [27], [28].  Other  work  using  the  p-version  finite 
elements  has  been  done  by  Babuska,  Craig,  Mandel,  Pitkaranta  [1]  and  Babuska, 
Griebel,  and  Pitkaranta  [2]. 

Much  work  using  domain  decomposition  has  focused  on  the  scalar  elliptic  problem, 
including  BJ0rstad  and  Widlund  [5], [6],  Bramble,  Pasciak,  and  Schatz  [7], [9], [10],  Chan 
and  Resasco  [11], [12],  Dryja  [16], [17],  Dryja  and  Widlund  [18], [19],  and  Widlund  [35]. 
The  work  by  Matsokin  and  Nepomnyaschikh  [29]  discusses  a  Schwarz  alternating 
algorithm  which  has  some  similarities  to  that  presented  in  this  paper. 

This  work  was  begun  while  visiting  the  University  of  Bergen,  Bergen,  Norway. 
I  would  like  to  thank  Petter  Bj0rstad  and  Anders  Hvidsten  for  their  assistance  in 
using  the  SESAM  code  and  for  many  useful  discussions.  I  would  also  like  to  thank 
VERITAS  SESAM  SYSTEMS,  H0vik,  Norway,  for  the  use  and  access  to  SESAM  for 
the  generation  of  test  stiffness  matrices  which  were  used  in  the  numerical  experiments 
reported  in  this  paper. 

2.  The  Discrete  Problem.  We  consider  a  second  order,  self  adjoint,  uniformly 
elliptic,  bilinear  form  an(u,r)  on  a  bounded  polyhedral  (polygonal)  domain  Q,  and, 
for  simplicity,  impose  a  homogeneous  Dirichlet  condition  on  dQ; 

au{u,v)  =  f{v),      yv  e  {H',(n))\      u  e  {H',iQ)y. 

Our  specific  applications  are  various  linear  elasticity  models,  where  q  is  from  2  to  6, 
cf.  [14].  Our  results  hold  also  for  the  scalar  case  q  —  1. 

We  perform  two  levels  of  triangulations  into  substructures  Q,  and  elements,  and 
assume  shape  regularity  and  that  the  substructures  and  elements  satisfy  the  usual 
rules  of  finite  element  triangulations;  see  e.g.  Ciarlet  [13].  V"{n)  C  {H^iQ.))"  and 
V''(fi)  C  {Hq{Q,))''  are  the  spaces  of  continuous,  piecewise  hnear  functions,  on  the 
two  triangulations,  which  vanish  on  the  boundary  dQ.  We  will  be  working  only  with 
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piecewise  lineax  finite  elements  but  note  that  the  theory  holds  also  for  higher  order 
elements. 

The  discrete  variational  problem  is  then  of  the  form 

(1)  anK,r,)  =  /K),    yv,ev\^),    u.eV^in). 

This  finite  dimensional  variational  problem  is  turned  into  a  linear  system  of  equations 
by  introducing  a  basis  {4>i}  for  the  space  V^.  We  use  the  standard  nodal  basis  functions 
in  V'*.  If  we  express  the  solution  as  Uh  =  ^Xi(}>i,  we  obtain  the  linear  system 

Kx  =  f. 

Here  x  is  the  vector  of  unknowns  x,, /  the  vector  of  components  f{(f>i),  and  A  the 
stiffness  matrix  A',j  =  aQ{(f>i^(f)j). 

Using  subscript  notations  explained  below,  for  any  block  matrix 

Kii     KiB 


Kb  I    Kbb 
with  A'//  nonsingular,  the  Schur  complement  of  K  is  given  by 

5  =  Kbb  -  KbiKjIKjb. 

This  is  the  block  which  remains  after  the  variables  associated  with  A'//  have  been 
eliminated  by  Gaussian  elimination.  We  give,  without  proof,  the  importcint  lemma 

Lemma  2.1.  If  K  is  -positive  definite  and  if  x  is  partitioned  in  the  same  manner 
as  K,x  =  {xj,xb)^,  then  for  any  given  xb, 

T  c  •       T  T' 

XgjXB  =  mm  a;   Ai. 

The  minimizing  extension  xj  is  referred  to  as  the  discrete  harmonic  extension.  This 
is  in  direct  2inaIogy  to  the  continuous  case  for  the  Laplacian.  The  discrete  harmonic 
extension  satisfies  the  relation  A//x/  +  KjbXb  =  0. 

With  each  substructure  Qkt  we  can  associate  the  portion  of  the  stiffness  matrix 
A'  arising  from  the  integration  over  that  substructure. 

Similarly,  the  subvector  of  unknowns  associated  with  this  substructure  is  called  x^'''. 
Since  we  are  working  with  nodal  basis  functions  the  support  of  </>,  associated  with  any 
node  in  the  interior  of  a  substructure  lies  completely  in  that  substructure.  Therefore 
there  axe  no  couplings  between  unknowns  associated  with  the  interiors  of  two  different 
substructures.  For  each  substructure,  we  partition  the  variables  into  those  associated 
with  a  node  in  the  interior  of  the  substructure,  x^  ',  and  those  on  the  boundary  of 
the  substructiure,  x^  .  To  simplify  notations,  we  pad  the  subvectors  x^*^^  and  subma- 
trices  A''*"'  by  zeros  so  that  the  following  equation  makes  sense.  This  is  called  the 
subassembly  process;  cf.   [19]. 


The  Schur  complement  of  the  global  stiffness  matrix  K  is  then  given  by 

=    E.5W4'- 

Each  of  the  Schur  complements,  5'*',  can  be  formed  independently  and  in  parallel. 
We  axe  now  left  with  the  still  large  linear  system 

Sxb  =  g. 

This  system  will  be  solved  using  preconditioned  conjugate  gradients.  The  precon- 
ditioner  is  constructed  using  the  ideas  developed  in  the  theory  of  additive  Schwarz 
methods. 

3.  Schwarz  Methods:  An  Abstract  Framework.  Suppose  that  we  wish  to 
solve  the  following  finite  dimensional  variational  problem 

(2)  a{u,v)  =  f{vl  yv£V,      ueV. 

Let  Vi  be  a  set  of  subspaces  of  V  so  that  V  =  Vi-\ 1-  Vyv  •  With  each  subspace  there 

is  a  corresponding  projection  operator  Pi,  which  is  the  projection  in  the  a{u,v)  inner 
product  onto  the  subspace  1^,  i.e. 

(3)  a{PiU,w)  =  a{u,w)^  f{w),  WweVi,       P,u  G  K-. 

PiU  can  be  determined  by  introducing  a  basis  {ipj  }  for  Vj  and  expanding  PiU  in  that 
basis,  PiU  =  J2j  ctj  4'j   ■  This  results  in  the  linear  system 

where  A'];   =  a(V'j   ,  V";    )?  and  /^'^  is  the  vector  defined  by  f{ipj). 

The  additive  Schwarz  method,  see  Dryja  and  Widlund  [19],  of  solving  equation  (2) 
is  to  introduce  an  auxiliary  problem 

(4)  Pu  =  Y,P,u  =  f, 

i 

which  has  the  same  solution  as  equation  (2).  Since  P,u  can  be  found  by  equation  (3) 
without  knowing  the  solution  of  (2),  we  first  compute  /  and  then  solve  equation  (4) 
using  the  conjugate  gradient  method. 

The  conjugate  gradient  method  can  be  an  effective  method  for  the  solution  of  a 
symmetric  (in  any  inner  product)  positive  definite  lineax  system.  The  decrease  in  the 
energy  norm  of  the  error  after  m  steps  can  be  bounded  by 


Jk{P)  -  1 
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where  k(P)  is  the  condition  number  of  the  matrix  P,  see  e.g.  Golub  and  Van  Loan 
[20]. 

The  reason  for  going  from  problem  (2)  to  problem  (4)  is  that  by  a  suitable  choice 
of  the  subspaces  Vi  we  cein  turn  a  lairge  ill-conditioned  system  into  a  very  well  con- 
ditioned problem  at  the  expense  of  solving  many  small  independent  linear  systems. 
The  following  two  lemmas  allow  us  to  develop  bounds  on  the  condition  number  of  P. 

Lemma  3.1.  Consider  the  undirected  graph  with  a  node  for  each  subspace  Vi,  and 
an  edge  between  node  i  and  node  j  iffViflVj  ^  0.  Let  p  be  the  number  of  colors  needed 
to  color  the  graph  so  that  no  two  nodes  connected  by  an  edge  have  the  same  color. 
Then 

Amax(P)   <  P- 


Proof  All  the  subspaces  for  a  pairticular  color  axe  disjoint,  hence  their  correspond- 
ing projection  operators  are  mutually  orthogonal.  Therefore  the  sum  of  the  projection 
operators  of  a  particular  color  is  itself  a  projection  operator.  P  then  consists  of  p  of 
these  composite  projection  operators  each  with  a  maximum  eigenvalue  of  one;  thus 

An^x(P)  <P.  ■ 

Lemma  3.2.  (Lions  [18], [24])  Assume  that  for  all  u  E  V,  there  exists  a  represen- 
tation u  =  Y^i  Ui  with  u,  e  Vi  such  that 

^a(ui,u,)  <  Coa(u,u) 
then 

Throughout  this  paper  c,  C  and  Co  represent  generic  constants  independent  of  h  and 
H. 


Proof. 


\u\l  =  ^a(u,u,)  =  ^a{u,P,Ui)  =  '^a{P,u,Ui). 


Therefore, 


w\i<ij:\pM.y^\E\^>\i) 


.|2^1/2 


By  the  assumption  of  the  lemma  and  a  property  of  projections, 

W\l  <  Co'E  l^."la  =  C^j:aiP,u,u)  =  CoyPu,u). 


4.   The  New  Algorithm.  The  variational  problem  we  are  solving  is 

ylSxB  =  ylg,         yyB  e  V\     xb  g  V". 

We  will  be  working  in  the  trace  space  {H^I'^iV))''  [21], [23].  The  space  V^  =  V^\r  is  a 
subspace  of  the  trace  space  {H^^'^(T))'' ,  and  the  quadratic  form  ygSws  is  a  discrete 
realization  of  the  (•,  ^//^/^(r)  inner  product  on  F  =  U5fi,.  This  can  easily  be  seen  by 
comparing  the  two  relationships, 

x^Sxb  =  niinx   Kx      and      (u,u);^i/2(r)  =  min(u,u)//i(n). 
^i  u|r=u 

For  our  algorithm  we  use  the  following  subspaces:  a  global  coarse  space,  and  spaces 
associated  with  overlapping  regions  of  F.  We  can  regard  the  boundaries  of  the  sub- 
structures as  consisting  of  three  pieces:  the  substructure  vertices,  the  edges  between 
substructure  vertices,  and  the  faces  of  the  substructures.  In  two  dimensions  there  are 
no  faces,  only  vertices  Eind  edges.  To  obtain  the  overlapping  regions  for  the  Schwarz 
method  we  introduce  the  faces  F^' .  We  choose  F-^-"  to  be  regions  consisting  of  an  edge 
and  an  overlap  of  order  H  onto  the  adjacent  faces.  The  F^*  axe  regions  consisting  of 
a  vertex  and  an  overlap  of  order  H  onto  all  adjacent  faces  and  edges.  We  constrain 
the  overlap  so  that  no  portion  of  F  is  covered  more  than  four  times. 
The  subspaces  of  V^  are  given  by, 

v;j=y''n(ifoT(r^"))% 

and 

{H]l'^{t))\t  C  F,  is  a  subspace  of  {H^I'^{V))'^ .  Its  norm  can  be  defined  by  \u\^m2^^^  = 

l"lHi/2(r)  where  u  \s  u  extended  by  zero  onto  F  \  F,  [21].  It  is  easy  to  see  that  each  of 
these  spaces  is  a  subspace  of  V^  and  that 

V'  =  (E  yk )  +  (E  K )  +  (l^  ^n)- 

The  projection  onto  each  'face'  subspace  is  given  by 

Pf,  =  Rf.^f,  Rf.S, 

where  Rf,  is  the  restriction  operator  which  returns  only  those  unknowns  which  are 
associated  with  F^'  and  Sf,  is  the  principal  minor  of  the  Schur  complement  S  which  is 
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associated  with  that  same  set  of  unknowns.  The  projections  onto  'edge'  and  'vertex' 
subspaces  are  formed  in  the  same  manner.  Naturally  the  restriction  operators  R 
need  never  be  exphcitly  represented,  instead  we  would  use  something  hke  hardware 
scatter-gather. 

For  the  coarse  space  F^,  the  projection  operator  is 

Ph  =  RIKh'RhS. 

The  operator  RJj  represents  linear  interpolation  from  V"  to  V^,  and  A'//  is  the  stiffness 
matrix  for  the  original  problem  treating  the  substructures  as  elements.  A'//  can  also 
be  obtained  by  making  a  partial  change  to  a  hierarchical  basis  of  each  S''''  so  that 
the  nodal  functions  associated  with  the  vertices  of  the  substructure  are  the  nodal 
functions  in  the  V"  space;  cf.  Smith  and  Widlund  [30].  The  principal  minor  of 
the  Schur  complement  which  is  associated  with  these  vertex  nodes  then  forms  the 
contribution  to  the  stiffness  matrix  A'//  from  the  given  substructure. 

If  we  exclude  from  the  algorithm  all  of  the  overlap  and  the  'vertex'  spaces  we  ob- 
tain the  iterative  substructuring  algorithm  presented  in  Dryja  and  Widlund  [19].  The 
iterative  substructuring  algorithm  is  also  very  similar  to  that  introduced  by  Bramble, 
Pasciak,  and  Schatz  [7].  In  two  dimensions  these  algorithms  have  condition  numbers 
which  grow  hke  {1  +  \og{H / h)y .  In  three  dimensions  the  condition  numbers  for  these 
algorithms  grow  faster  than  (H/h),  [8], [16]. 

We  now  present  a  detailed  description  of  the  algorithm.  We  note  that  many 
opportunities  exist  for  parallelism  between  and  within  each  step. 

Iterative  substructuring/ additive  Schwarz  Algorithm 

1.  Form  the  stiffness  matrices 

/    L-(j)       tAj)  \ 

\  ^IB        ^BB  J 

for  each  substructure  by  integration. 

2.  Factor  Kj]  for  each  substructure. 

3.  Form  the  Schur  complements  5^^)  =  K^Jl  -  A"/'/a'}?"'A'}2. 

4.  Form  ajid  factor  the  coEirse  stiffness  matrix  A'//. 

5.  Form,  by  subassembly,  SpiiSsi,  and  5v;. 

6.  Factor  Sf,,Se,i  and  Sv, ■ 

7.  Form  the  right  hand  sides  b^^^  for  each  substructure  by  integration. 

8.  Modify  the  right  hand  sides;  fc^'  =  6^^  -  K^'^'' K\i'>~' b\'\ 

9.  Solve  Sxb  =  ^B  using  a  preconditioned  conjugate  gradient  method  with  the 
preconditioner 

S~^  =  RffK^^RH  +  Yl  Rf,^f!Rf.  +  ^  Re,Se]Re,  +  XI  ^v^^vl^v^- 

J  k 

10.   Form  the  right  hand  sides  for  the  problems  on  the  interiors  of  the  subdomains 

i,(i)_/,0)_;^(i)-'E'0)3.(i) 

Oj       —    Oj       —  ri.jj         I\.jgXg    . 


11.  Solve  for  the  interior  unknowns  x)    =  Kj]     b/' . 
Steps  1  though  6  can  be  regarded  as  a  preprocessing  stage  independent  of  the  partic- 
ulcir  loads;  they  need  not  be  repeated  for  different  loads. 

The  main  result  of  this  paper  is  now  given  in 

Theorem  4.1.  The  condition  number  of  the  preconditioned  linear  system  is 
bounded  independently  of  the  size  of  the  substructures  H  and  the  size  of  the  elements 
h,  i.e. 

k{P)  =  k(5-^5)  <  C. 

Proof  We  note  that  no  node  on  T  is  contained  in  more  than  4  of  the  regions 
pF,  ^  pBj  ^j^^  pi4  Therefore  no  element  of  the  space  V^  can  belong  to  more  than  5 
subspaces  (the  four  above  and  V^);  hence  by  Lemma  3.1  we  have  Ainax(-P)  ^  5. 

We  first  use  the  uniform  ellipticity  of  our  bilinear  form  an(u,  ti), 

c|"l//i(n)  <  an(u,u)  <  C|u|^:(n), 

to  meike  it  possible  to  work  in  the  (if^(n))'  and  (if^/^(r))'  norms  instead  of  the 
equivalent  norms  induced  by  an(u,  u).  To  obtain  a  lower  bound,  we  must  demonstrate 
that  for  all  u'^  £  V^,  there  exists  a  representation 

t         j         k 

so  that 

\u    |//i/2(r)  +  ^«  I^JE,  lH'/2(r)  "*"  ^j  I^F,  l//i/2(r)  "*"  ^*:  1^ v)i '//^/^(r) 
(5) 

—  Co\u  |//i/2(r)» 

where  Co  is  independent  oi  u^,h  and  H. 

We  construct  this  representation  as  follows.  We  extend  u  as  a  discrete  har- 
monic function  into  Q.  ,  with  u^\r  =  u^.  In  addition  we  extend  the  boundary  regions 
pE,  ^pF,  pVfc  jj^^Q  ^YiQ  neighboring  substructures.  This  results  in  a  large  collection  of 
overlapping  regions.  Now  apply  the  result  of  Appendix  A  to  each  component  of  the 
vector  valued  function  u'^  separately  to  obtain  a  representation  of  u'', 


u''  =  ^.«  +  ^4.  +  ^4^  +  ^<, 


which  satisfies 

t  j  k 

We  now  define  our  representation  by 

u"  =  u"\r,      4.  =  4jr, 


The  use  of  the  trace  theorem  [21], 

c|u|f/i/2(r)  <  \u\hhq) 

and  an  extension  theorem  [34], 

\u''\hhq)  <  C\u  \HU2^^) 

then  gives  the  bound  (5).  We  then  apply  Lemma  (3.2)  to  conclude  the  proof  of  the 
theorem.  ■ 

5.  Numerical  Results.  We  have  performed  numerous  experiments  with  prob- 
lems in  two  dimensions.  The  problems  considered  are 

•  The  Laplacian  using  the  usual  five  point  stencil; 

•  The  equations  of  linear  elasticity  using  4  node  square  membrane  elements 
with  2  degrees  of  freedom  per  node;  cf.  [3]; 

•  The  equations  of  linear  elasticity  using  4  node  square  shell  elements  with  3 
degrees  of  freedom  per  node;  cf.  [3]. 

Experiments  have  been  performed  on  square  and  L-shaped  regions;  since  there  was 
no  appreciable  difference  between  the  two  cases  results  are  only  given  for  the  square 
regions.  The  substructures  are  squares.  For  the  elasticity  problems  the  stiffness  matri- 
ces were  generated  using  the  SESAM  code  [3],  a  large,  reliable  commercial  structural 
analysis  code,  using  a  Poisson  ratio  of  .3. 

The  experiments  were  run  twice,  once  using  all  the  subspax:es  as  indicated  in  the 
algorithm  and  once  excluding  the  'vertex'  spaces.  As  expected  the  condition  number 
remains  bounded  by  a  constant  independent  of  H  and  h  when  all  the  spaces  were 
included.  When  the  'vertex'  spaces  were  excluded  the  condition  number  appeairs  to 
grow  like  (1  +  \og{H / h)y ,  as  expected. 

The  selection  of  an  appropriate  stopping  condition  for  the  preconditioned  con- 
jugate gradient  method  is  crucizd.  A  stopping  criteria  based  only  on  a  norm  of  the 
residual  can  make  comparisons  between  preconditioned  and  unpreconditioned  results 
misleading  since  the  eigenvalues  of  the  two  operators  can  be  of  completely  differ- 
ent orders  of  magnitude.  For  instance,  for  elasticity  problems  the  eigenvalues  of  the 
original  stiiTness  matrices  can  be  of  order  10^^  while  the  eigenvalues  of  the  precondi- 
tioned problems  generally  are  of  order  1 .  We  have  therefore  chosen  to  use  the  stopping 
condition 

II      .J     111      /  e||approx.  solution||£,2 
||residual||/,2  <  ; . 

^imn{S~^ S)  is  calculated  using  the  Lanczos  method  at  very  little  extra  expense.  We 
have  chosen  to  use  e  =  10"^;  this  assures  that  roughly  five  digits  of  the  solution  are 
correct  and  not  many  more,  regardless  of  the  preconditioner  used. 


Number        Nodes     Number  of  No  Without 

of  along      Unknowns     Preconditioner       'Vertex' 

Subdomains      Edge  on  T  Spaces 


With 

'Vertex' 

Spaces 


16 


3 

7 
15 

31 


81 
177 
369 
753 


35.26 

75.10 

155 

315 


14 
24 
37 
51 


4.92  7 

7.59  9 

10.89  10 

14.84  11 


2.45  6 

2.55  7 

2.82  7 

2.99  7 


64 


3 

7 

15 

31 


385 
833 

1729 
3521 


137 
290 
598 
1216 


32 

49 
70 


5.35  9 

8.19  10 

11.54  12 

15.62  13 


2.60  8 

2.68  8 

2.78  8 

2.87  8 


256 


3 
7 

15 
31 


1665 
3585 
7425 
15105 


545 
1151 
2372 
4766 


62 

91 

* 


5.46  9 

8.27  10 

11.77  12 

15.90  13 


2.63  8 

2.70  7 

2.80  7 

2.89  7 


Table  1 
Condition  Numbers  and  Iteration  Counts  for  the  Laplacian 


Overlap  in  nodes 

0 

1 

2 

3 

4 

5 

6 

7 

8 

Condition  number 

15.62 

4.49 

4.01 

3.78 

3.52 

3.38 

3.01 

2.92 

2.81 

Iterations 

13 

8 

8 

8 

8 

8 

8 

8 

8 

Table  2 
Condition  Number  as  Function  of  Overlap  for  the  Laplacian 


No  Preconditioning 


Iter. 


t 

1 
2 
3 

4 
5 
6 

7 

8 

9 

10 

11 

12 


Without  'Vertex' 
Spaces 


With  'Vertex' 
Spaces 


lleib 

9.8  X  10-^ 
9.6  X  10-^ 

9.5  X  10-1 

9.3  X  10-1 
9.2  X  10-1 

9.1  X  IQ-i 

8.9  X  10-1 

8.8  X  10-1 

8.6  X  IQ-i 

8.4  X  10-1 

8.2  X  10-1 

7.9  X  10-1 


.99 
.98 
.98 
.98 
.98 
.98 
.98 
.98 
.98 
.98 
.98 
.98 


3.0  X  10-2 
9.9  X  10-3 
2.2  X  10-3 

1.2  X  10-3 

6.3  X  10-^ 

4.5  X  10-^ 
1.0  X  10-^ 

3.6  X  10-5 
9.5  X  10-^ 
6.3  X  10-^ 

7.7  X  10-^ 

3.0  X  10-^ 
Table  3 


( 


.03 
.10 
.13 
.19 
.23 
.28 
.27 
.28 
.28 
.30 
.34 
.35 


1/. 


3.0  X  10-2 

1.0  X  10-2 

4.1  X  10-3 
9.0  X  lO-'' 
7.0  X  10-5 
2.5  X  10-5 
9.4  X  10-^ 

3.3  X  10-^ 

3.8  X  10-^ 
6.7  X  10-« 

2.4  X  10-« 

5.9  X  10-^ 


.03 


.10 
.16 
.17 
.15 
.17 
.19 
.21 
.19 
.19 
.20 
.21 


Errors  and  Convergence  Rates  for  Laplacian 


Number 

Nodes 

Number  of 

No 

Without 

With 

of 

along 

Unknowns 

Preconditioner 

'Vertex' 

'Vertex' 

Subdomains 

Edge 

on  r 

Spaces 

Spaces 

16 

3 

162 

22.16 

19 

10.52 

15 

3.51     10 

7 

354 

46.63 

28 

14.84 

17 

3.51     10 

15 

738 

96.34 

41 

19.83 

19 

3.56     10 

31 

1506 

196 

60 

25.51 

20 

3.62     10 

64 

3 

770 

84.82 

37 

12.13 

18 

3.85     10 

7 

1666 

178 

55 

16.92 

19 

3.85     10 

15 

3458 

368 

79 

22.37 

22 

3.84     10 

31 

7042 

747 

* 

28.52 

25 

3.89     10 

256 

3 

3330 

334 

75 

12.42 

18 

3.91     10 

7 

7170 

705 

* 

17.31 

19 

3.90     10 

15 

14850 

1453 

+ 

22.88 

22 

3.89     10 

31 

30210 

2921 

* 

29.15 

25 

3.94     10 

Table  4 
Condition  Numbers  and  Iteration  Counts  for  Membrane  Elements 


Overlap  in  nodes 

0 

1 

2 

3 

4 

5 

6 

7 

8 

Condition  number 

28.55 

5.59 

4.85 

4.36 

4.01 

3.89 

3.88 

3.89 

3.89 

Iterations 

25 

12 

11 

10 

10 

10 

10 

10 

10 

Table  5 
Condition  Number  as  Function  of  Overlap  for  Membrane  Elements 


Iter. 

No  Preconditioning 

Without  'Vertex' 
Spaces 

With  'Vertex' 
Spaces 

i 

1 

9.7  X  10-^ 

(111 
M|e' 

.98 

1  e-  \l^ 
8.0  X  10-2 

(ih 

•11,2X1/. 

.08 

8.7  X  10-2 

Mk 

'11,2X1/,- 
"1112'' 

.09 

2 

9.5  X  10-^ 

.98 

3.7  X  10-2 

.19 

4.2  X  10-2 

.20 

3 

9.3  X  10-1 

.98 

1.1  X  10-2 

.22 

1.6  X  10-2 

.25 

4 

9.2  X  10-1 

.98 

6.3  X  10-3 

.28 

3.9  X  10-3 

.25 

5 

9.0  X  10-1 

.98 

3.7  X  10-3 

.33 

1.1  X  10-3 

.26 

6 

8.8  X  IQ-i 

.98 

3.4  X  10-3 

.39 

5.1  X  10-1 

.28 

7 

8.6  X  IQ-i 

.98 

2.0  X  10-3 

.41 

1.6  X  10-1 

.29 

8 

8.3  X  10-1 

.98 

1.4  X  10-3 

.44 

4.4  X  10-^ 

.29 

9 

8.1  X  10-1 

.98 

5.1  X  10-1 

.43 

1.7  X  10-^ 

.30 

10 

7.8  X  IQ-i 

.98 

2.3  X  10-" 

.43 

5.4  X  10-^ 

.30 

11 

7.5  X  10-1 

.97 

2.5  X  10-" 

.47 

1.5  X  10-^ 

.30 

12 

7.2  X  10-1 

Errors 

.97 

and  Con 

2.6  X  10-1 
Table  6 
vergence  Rates  j 

.50 

for  Membrane 

5.1  X  10-' 

Elements 

.30 

10 


Number 

Nodes 

Number  of 

No 

Without 

With 

of 

along 

Unknowns 

Preconditioner 

'Vertex' 

'Vertex' 

Subdomains 

Edge 

on  r 

Spaces 

Spaces 

16 

3 

243 

452 

47 

10.51 

17 

3.48     10 

7 

531 

442 

66 

14.83 

18 

3.48     10 

15 

1107 

970 

* 

19.81 

21 

3.53     10 

31 

2259 

* 

* 

25.48 

23 

3.60     10 

64 

3 

1155 

1751 

* 

11.99 

18 

3.72     10 

7 

2499 

1707 

* 

15.93 

19 

3.74     10 

15 

5187 

3800 

* 

22.12 

23 

3.83     10 

31 

10563 

* 

* 

28.24 

26 

3.88     10 

256 

3 

4995 

* 

* 

11.63 

17 

3.73     10 

7 

10755 

* 

♦ 

16.19 

19 

3.81     10 

15 

22275 

* 

+ 

22.58 

22 

3.85     10 

Table  7 
Condition  Numbers  and  Iteration  Counts  for  Shell  Elements 


Overlap  in  nodes 

0 

1 

2 

3 

4 

5 

6 

7 

8 

Condition  number 

28.24 

5.62 

4.81 

4.33 

3.98 

3.87 

3.87 

3.88 

3.88 

Iterations 

23 

12 

11 

10 

10 

10 

10 

10 

10 

Table  8 
Condition  Number  as  Function  of  Overlap  for  Shell  Elements 


No  Preconditioning 

Without 

'Vertex' 

With  '^ 

V^ert. 

5X' 

Iter. 

Spaces 

Spaces 

i 
1 

* 

* 

?" 

6.8  X  10-2 

'\\l2> 
.07 

7.4  X  10-2 

(^ 

.07 

2 

* 

* 

3.2  X  10-2 

.18 

3.6  X  10-2 

.19 

3 

* 

* 

9.4  X  10-3 

.21 

1.4  X  10-2 

.24 

4 

* 

* 

5.3  X  10-3 

.27 

3.4  X  10-3 

.24 

5 

* 

* 

3.1  X  10-3 

.32 

1.0  X  10-3 

.25 

6 

* 

* 

3.1  X  10-3 

.38 

4.7  X  10-" 

.28 

7 

* 

* 

1.9  X  10-3 

.41 

1.4  X  10-" 

.28 

8 

* 

* 

1.4  X  10-3 

.44 

3.9  X  10-^ 

.28 

9 

+ 

* 

5.5  X  10-" 

.43 

1.6  X  10-^ 

.29 

10 

* 

* 

3.3  X  lO-'' 

.45 

4.9  X  10-^ 

.29 

11 

* 

* 

2.5  X  10-" 

.47 

1.3  X  10-^ 

.29 

12 

* 

* 

2.5  X  10-" 

.50 

4.6  X  10-^ 

.30 

13 

* 

* 

2.4  X  10-" 

Table  9 

.53 

1.5  X  10-^ 

.30 
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In  Table  1  the  experiments  are  conducted  for  the  Laplacian.  The  overlap  of  the 
'vertex'  spaces  onto  the  'edge'  spaces  is  chosen  to  be  /f/4.  In  Table  2,  we  examine  the 
effect  of  varying  the  amount  of  overlap  of  the  'vertex'  spaces  onto  the  'edge'  spaces 
for  the  case  with  64  substructures  and  31  nodes  along  the  edge  of  each  substructure. 
We  see  that  the  overlap  is  very  important  but  a  small  overlap  has  almost  as  much 
effect  as  a  larger  overlap.  We  give  a  sample  of  the  convergence  behavior  in  Table  3, 
showing  the  discrete  L^  norm  of  the  error  as  a  function  of  the  number  of  iterations. 
This  is  again  for  the  case  of  64  substructures  and  31  nodes  along  the  edge  of  each 
substructure.  These  tables  are  repeated  for  the  linear  elasticity  problems. 

A.  Appendix:  A  Partitioning  Result.  For  two  or  three  dimensions  let  Q. 
be  a  polyhedral  (polygonal)  domain  which  has  been  triangulated  into  substructures 
which  are  shape  regulcir,  with  diameter  0(H).  Continue  the  triangulation  to  obtain  a 
triangulation  with  elements  of  diameter  0(h).  Furthermore  assume  that  Q  has  been 
covered  with  N  shape  regular  overlapping  regions  fi,,  (not  necessarily  related  to  the 
coarse  triangulation  above)  each  with  a  dicimeter  0(H).,  each  of  which  overlaps  all 
its  neighbors  with  an  overlap  of  0(H).  Let  V"(9.)  C  H^(n)  and  V''(n)  C  H^(n)  be 
the  spaces  of  continuous,  piecewise  linear  functions,  on  the  two  triangulations,  which 
vanish  on  the  boundary  dO,.  We  then  construct  the  following  spaces 

Vo^  =  V",      V,'' =  V'' n  H^(Q,). 

The  following  theorem  is  a  variation  of  a  result  given  in  Dryja  and  Widlund  [19]. 
Theorem  A.l.  For  all  u^  e  V"  there  exists  u-"  G  V;''  with  u''  =  XI^o  "I'  ^"c/i  that 

N 

where  Co  is  independent  of  u\h  and  H. 

Proof.  From  Strang  [32],  we  know  that  there  exists  a  linear  map  Ih  '•  V^  —>  V^ 
which  satisfies 

(6)  \\n'-iHu''\\l^(^)<CH'\u%.(^) 
and 

(7)  |u''-///"''lHMn)<C'|u''|^:(n). 
We  then  define 

h  ...h         ..h        ^,h 


and 


<  =  ///u^       zx;"  =  u'^  -  < 


uf  =  h(e,w''). 
12 


Ih  is  the  lineax  interpolation  operator  onto  the  space  V'^  and  the  6,  form  a  partition 
of  unity  with  6i  e  C^{9,i),0  <  6,  <  1  and  ^^j  ^,  =  1.  Since  h  is  a  Unear  operator,  it 
is  immediate  that 

N 

Because  of  the  generous  overlap  between  subregions,  we  ccin  insure  that  the  gradients 
of  6i  axe  well  behaved.  That  is,  6i  can  be  constructed  so  that  its  gradients  never 
grow  faster  than  |V^,||oo  <  C/H^.  If  we  let  K  represent  any  single  element  in  the 
triangulation  this  can  be  expressed  as 

(8)  \\o,-e,\\l^^j,)<c{h/HY. 

Here  6i  is  the  average  of  6,  on  element  K. 

We  now  estimate  the  H^  norm  of  u'l  over  a  single  element. 

which  can  be  bounded  using  an  inverse  inequality  by 

\^'l\mK)   <  m^^%^^K)  +  Ch-'\\h{e,-e,)w'^)\\l^^.y 

We  now  use  equation  (8)  and  the  trivial  inequality  ||^,||l°°  <  1  to  obtain 

Since  a  finite  bounded  number  of  u,^  axe  nonzero  for  any  element  K,  we  obtain,  when 
summing  over  i, 

t=l 
Next  sum  over  the  elements  K, 

To  finish  the  argument,  we  use  equations  (6)  and  (7)  to  obtain 

N 
t=0 
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